function [e] = fc_wk4(p);

global pressure flow T

% function 4-el WK-model
% used to estimate L and C via 4-el WK model
% method: put flow in m-el model
% minimise difference between measured and computed pressure
% p(1)=TVR
% p(2)=Z0
% p(3)=C
% p(4)=L


aantal=length(pressure);

fftQ=fft(flow);

R_WK4=p(1);
Z0_WK4=p(2);
C_WK4=p(3);
L_WK4=p(4);

% application 4element WK model: compute Pressure response of 3-el WK model to flow

% step 1: transform Matlab fft flow into harmonics

fftQ(1)=fftQ(1)/aantal;

for j=2:aantal/2
    fftQ(j)=2*fftQ(j)/aantal;
end

% step 2: compute input impedance 3-el WK model and harmonics P_WK4 

for j=1:aantal/2
    freq(j)=(j-1)/T;
    omega=2*pi*freq(j);
    Z_WK4(j)=i*omega*Z0_WK4*L_WK4/(Z0_WK4+i*omega*L_WK4)+R_WK4/(1+i*omega*R_WK4*C_WK4);
    fftP_WK4(j)=Z_WK4(j)*fftQ(j);
end

% step 3: compute matlab harmonics

fftP_WK4(1)=fftP_WK4(1)*aantal;

% check for even or uneven number

if rem(aantal,2)==0  
	for j=1:aantal/2-1
	    fftP_WK4(j+1)=fftP_WK4(j+1)*aantal/2;
	    fftP_WK4(aantal-j+1)=conj(fftP_WK4(j+1));
	end;
	fftP_WK4(aantal/2)=fftP_WK4(j+1)*aantal/2;
else
	for j=1:ceil(aantal/2)-1
	    fftP_WK4(j+1)=fftP_WK4(j+1)*aantal/2;
	    fftP_WK4(aantal-j+1)=conj(fftP_WK4(j+1));
	end;
end;


% step 4: compute inverse matlab transform

P_WK4=real(ifft(fftP_WK4));

if size(pressure)==size(P_WK4)
   verschil=(pressure-P_WK4).*(pressure-P_WK4);
else
   verschil=(pressure-P_WK4').*(pressure-P_WK4');
end;   

e=sum(verschil);



